function [A,B,Pi] = rndinit(n_s,n_o,uA,uB,uPi,N,total_len)

wa = zeros(n_s, n_s);
wb = zeros(n_s, n_o);

wpi = dirrnd(uPi,1)*N;
for k=1:n_s, % loop over hidden states
    wa(k,:) = dirrnd(uA(k,:),1)*total_len;
    wb(k,:) = dirrnd(uB(k,:),1)*total_len;
end;
Wa = wa + uA; % posterior is data counts plus prior.
Wb = wb + uB;
Wpi = wpi + uPi;

A = Wa./  repmat( sum(Wa,2) ,[1 n_s])  ;
B = Wb./  repmat( sum(Wb,2) ,[1 n_o])  ;
Pi = Wpi/  sum(Wpi,2)   ;

end

